####Model results####
##a) policy (table X29, models A99-A102)
i4a_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_policy_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
i4a_m1g_t1_sc_cse <- data.frame(cluster.se(i4a_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
i4a_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_policy_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
i4a_m1g_t1_sf_cse <- data.frame(cluster.se(i4a_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
i6a_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_policy_territory_t","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6a_m1d_t1_sb_cse <- data.frame(cluster.se(i6a_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
i6a_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_policy_territory_t*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6a_m2d_t1_sb_cse <- data.frame(cluster.se(i6a_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(i4a_m1g_t1_sc, i4a_m1g_t1_sf, i6a_m1d_t1_sb, i6a_m2d_t1_sb, se=c(i4a_m1g_t1_sc_cse, i4a_m1g_t1_sf_cse, i6a_m1d_t1_sb_cse, i6a_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A99)", "Min. (model A100)","Maj./min. dyad (model A101)","Maj./min. dyad (model A102)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X29. Additional results: policy autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (policy)","Territorial autonomy (policy) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex29.txt")
stargazer(i4a_m1g_t1_sc, i4a_m1g_t1_sf, i6a_m1d_t1_sb, i6a_m2d_t1_sb, se=c(i4a_m1g_t1_sc_cse, i4a_m1g_t1_sf_cse, i6a_m1d_t1_sb_cse, i6a_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A99)", "Min. (model A100)","Maj./min. dyad (model A101)","Maj./min. dyad (model A102)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X29. Additional results: policy autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (policy)","Territorial autonomy (policy) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex29.html")
##b) fiscal (table X30, models A103-A106)
i4b_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_fiscal_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
i4b_m1g_t1_sc_cse <- data.frame(cluster.se(i4b_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
i4b_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_fiscal_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
i4b_m1g_t1_sf_cse <- data.frame(cluster.se(i4b_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
i6b_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_fiscal_territory_t","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6b_m1d_t1_sb_cse <- data.frame(cluster.se(i6b_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
i6b_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_fiscal_territory_t*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6b_m2d_t1_sb_cse <- data.frame(cluster.se(i6b_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(i4b_m1g_t1_sc, i4b_m1g_t1_sf, i6b_m1d_t1_sb, i6b_m2d_t1_sb, se=c(i4b_m1g_t1_sc_cse, i4b_m1g_t1_sf_cse, i6b_m1d_t1_sb_cse, i6b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A103)", "Min. (model A104)","Maj./min. dyad (model A105)","Maj./min. dyad (model A106)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X30. Additional results: fiscal autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (fiscal)","Territorial autonomy (fiscal) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex30.txt")
stargazer(i4b_m1g_t1_sc, i4b_m1g_t1_sf, i6b_m1d_t1_sb, i6b_m2d_t1_sb, se=c(i4b_m1g_t1_sc_cse, i4b_m1g_t1_sf_cse, i6b_m1d_t1_sb_cse, i6b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A103)", "Min. (model A104)","Maj./min. dyad (model A105)","Maj./min. dyad (model A106)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X30. Additional results: fiscal autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (fiscal)","Territorial autonomy (fiscal) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex30.html")
##c) political (table X31, models A107-A110)
i4c_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_political_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
i4c_m1g_t1_sc_cse <- data.frame(cluster.se(i4c_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
i4c_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_political_territory_t","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
i4c_m1g_t1_sf_cse <- data.frame(cluster.se(i4c_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
i6c_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_political_territory_t","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6c_m1d_t1_sb_cse <- data.frame(cluster.se(i6c_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
i6c_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_political_territory_t*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
i6c_m2d_t1_sb_cse <- data.frame(cluster.se(i6c_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(i4c_m1g_t1_sc, i4c_m1g_t1_sf, i6c_m1d_t1_sb, i6c_m2d_t1_sb, se=c(i4c_m1g_t1_sc_cse, i4c_m1g_t1_sf_cse, i6c_m1d_t1_sb_cse, i6c_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A107)", "Min. (model A108)","Maj./min. dyad (model A109)","Maj./min. dyad (model A110)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X31. Additional results: political autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (political)","Territorial autonomy (political) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex31.txt")
stargazer(i4c_m1g_t1_sc, i4c_m1g_t1_sf, i6c_m1d_t1_sb, i6c_m2d_t1_sb, se=c(i4c_m1g_t1_sc_cse, i4c_m1g_t1_sf_cse, i6c_m1d_t1_sb_cse, i6c_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A107)", "Min. (model A108)","Maj./min. dyad (model A109)","Maj./min. dyad (model A110)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X31. Additional results: political autonomy component only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (political)","Territorial autonomy (political) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex31.html")

####Figure A17####
##a) policy (civil)
me_i4a_m1g_t1_sc <- margins_summary(i4a_m1g_t1_sc, variables = "sa_policy_territory_t", vcov = cluster.vcov(i4a_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_i4a_m1g_t1_sc$term <- "second-order maj."
me_i4a_m1g_t1_sf <- margins_summary(i4a_m1g_t1_sf, variables = "sa_policy_territory_t", vcov = cluster.vcov(i4a_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_i4a_m1g_t1_sf$term <- "second-order min."
me_i4a_m1g <- rbind.fill(me_i4a_m1g_t1_sc, me_i4a_m1g_t1_sf)
me_i4a_m1g$term <- as.factor(me_i4a_m1g$term)
me_i4a_m1g$model <- "policy autonomy"
me_i4a_m1g$estimate <- me_i4a_m1g$AME
me_i4a_m1g$conf.low <- me_i4a_m1g$lower
me_i4a_m1g$conf.high <- me_i4a_m1g$upper
gc()
##b) policy (communal)
me_i6a_m1d_t1_sb <- margins_summary(i6a_m1d_t1_sb, variables = "sa_policy_territory_t", vcov = cluster.vcov(i6a_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6a_m1d_t1_sb$term <- "all"
me_i6a_m1d_t1_sb$model <- "policy autonomy"
me_i6a_m2d_t1_sb <- margins_summary(i6a_m2d_t1_sb, variables = "sa_policy_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(i6a_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6a_m2d_t1_sb$term <- ifelse(me_i6a_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_i6a_m2d_t1_sb$model <- "policy autonomy"
me_i6a_m13d <- rbind.fill(me_i6a_m1d_t1_sb, me_i6a_m2d_t1_sb)
me_i6a_m13d <- subset(me_i6a_m13d, term != "other")
me_i6a_m13d$term <- as.factor(me_i6a_m13d$term)
me_i6a_m13d$estimate <- me_i6a_m13d$AME
me_i6a_m13d$conf.low <- me_i6a_m13d$lower
me_i6a_m13d$conf.high <- me_i6a_m13d$upper
##c) fiscal (civil)
me_i4b_m1g_t1_sc <- margins_summary(i4b_m1g_t1_sc, variables = "sa_fiscal_territory_t", vcov = cluster.vcov(i4b_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_i4b_m1g_t1_sc$term <- "second-order maj."
me_i4b_m1g_t1_sf <- margins_summary(i4b_m1g_t1_sf, variables = "sa_fiscal_territory_t", vcov = cluster.vcov(i4b_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_i4b_m1g_t1_sf$term <- "second-order min."
me_i4b_m1g <- rbind.fill(me_i4b_m1g_t1_sc, me_i4b_m1g_t1_sf)
me_i4b_m1g$term <- as.factor(me_i4b_m1g$term)
me_i4b_m1g$model <- "fiscal autonomy"
me_i4b_m1g$estimate <- me_i4b_m1g$AME
me_i4b_m1g$conf.low <- me_i4b_m1g$lower
me_i4b_m1g$conf.high <- me_i4b_m1g$upper
##d) fiscal (communal)
me_i6b_m1d_t1_sb <- margins_summary(i6b_m1d_t1_sb, variables = "sa_fiscal_territory_t", vcov = cluster.vcov(i6b_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6b_m1d_t1_sb$term <- "all"
me_i6b_m1d_t1_sb$model <- "fiscal autonomy"
me_i6b_m2d_t1_sb <- margins_summary(i6b_m2d_t1_sb, variables = "sa_fiscal_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(i6b_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6b_m2d_t1_sb$term <- ifelse(me_i6b_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_i6b_m2d_t1_sb$model <- "fiscal autonomy"
me_i6b_m13d <- rbind.fill(me_i6b_m1d_t1_sb, me_i6b_m2d_t1_sb)
me_i6b_m13d <- subset(me_i6b_m13d, term != "other")
me_i6b_m13d$term <- as.factor(me_i6b_m13d$term)
me_i6b_m13d$estimate <- me_i6b_m13d$AME
me_i6b_m13d$conf.low <- me_i6b_m13d$lower
me_i6b_m13d$conf.high <- me_i6b_m13d$upper
##e) political (civil)
me_i4c_m1g_t1_sc <- margins_summary(i4c_m1g_t1_sc, variables = "sa_political_territory_t", vcov = cluster.vcov(i4c_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_i4c_m1g_t1_sc$term <- "second-order maj."
me_i4c_m1g_t1_sf <- margins_summary(i4c_m1g_t1_sf, variables = "sa_political_territory_t", vcov = cluster.vcov(i4c_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_i4c_m1g_t1_sf$term <- "second-order min."
me_i4c_m1g <- rbind.fill(me_i4c_m1g_t1_sc, me_i4c_m1g_t1_sf)
me_i4c_m1g$term <- as.factor(me_i4c_m1g$term)
me_i4c_m1g$model <- "political autonomy"
me_i4c_m1g$estimate <- me_i4c_m1g$AME
me_i4c_m1g$conf.low <- me_i4c_m1g$lower
me_i4c_m1g$conf.high <- me_i4c_m1g$upper
##f) political (communal)
me_i6c_m1d_t1_sb <- margins_summary(i6c_m1d_t1_sb, variables = "sa_political_territory_t", vcov = cluster.vcov(i6c_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6c_m1d_t1_sb$term <- "all"
me_i6c_m1d_t1_sb$model <- "political autonomy"
me_i6c_m2d_t1_sb <- margins_summary(i6c_m2d_t1_sb, variables = "sa_political_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(i6c_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_i6c_m2d_t1_sb$term <- ifelse(me_i6c_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_i6c_m2d_t1_sb$model <- "political autonomy"
me_i6c_m13d <- rbind.fill(me_i6c_m1d_t1_sb, me_i6c_m2d_t1_sb)
me_i6c_m13d <- subset(me_i6c_m13d, term != "other")
me_i6c_m13d$term <- as.factor(me_i6c_m13d$term)
me_i6c_m13d$estimate <- me_i6c_m13d$AME
me_i6c_m13d$conf.low <- me_i6c_m13d$lower
me_i6c_m13d$conf.high <- me_i6c_m13d$upper
##g) combined (civil)
figurea17_plot_data_g <- rbind.fill(me_i4a_m1g,me_i4b_m1g, me_i4c_m1g)
figurea17_plot_data_g$model <- as.factor(figurea17_plot_data_g$model)
figurea17_plot_g <- dwplot(figurea17_plot_data_g,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence") +
  scale_color_manual(name = "Coefficient for:", values = c("#1b9e77","#d95f02","#7570b3"))+ scale_linetype_manual(name = "Coefficient for:", values = c(3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+ #+
  guides(shape = guide_legend(nrow=1,"model",reverse=TRUE), colour = guide_legend(nrow=1,"model",reverse=TRUE),linetype = guide_legend(nrow=1,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##h) combined (communal)
figurea17_plot_data_d <- rbind.fill(me_i6a_m13d,me_i6b_m13d,me_i6c_m13d)
figurea17_plot_data_d$model <- as.factor(figurea17_plot_data_d$model)
figurea17_plot_d <- dwplot(figurea17_plot_data_d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) communal violence maj./min. dyad") +
  scale_color_manual(name = "Coefficient for:", values = c("#1b9e77","#d95f02","#7570b3"))+ 
  scale_linetype_manual(name = "Coefficient for:", values = c(3,2,1))+ 
  scale_shape_manual(name = "Coefficient for:", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type")+ #+
  guides(shape = guide_legend(nrow=1,"model",reverse=TRUE), colour = guide_legend(nrow=1,"model",reverse=TRUE),linetype = guide_legend(nrow=1,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##i) combined
figurea17 <- grid_arrange_shared_legend(figurea17_plot_g, figurea17_plot_d, ncol=2, nrow=1)
ggsave(figurea17, file='../figures/figurea17.pdf', width = 14, height = 5, units="cm",dpi=1000)